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Abstract — We develop mask iterative hard thresholding al- 
gorithms (mask IHT and mask DORE) for sparse image re- 
construction of objects with known contour. The measurements 
follow a noisy underdetermined linear model common in the 
compressive sampling literature. Assuming that the contour of 
the object that we wish to reconstruct is known and that the 
signal outside the contour is zero, we formulate a constrained 
residual squared error minimization problem that incorporates 
both the geometric information (i.e. the knowledge of the object's 
contour) and the signal sparsity constraint. We first introduce a 
mask IHT method that aims at solving this minimization problem 
and guarantees monotonically non-increasing residual squared 
error for a given signal sparsity level. We then propose a double 
overrelaxation scheme for accelerating the convergence of the 
mask IHT algorithm. We also apply convex mask reconstruction 
approaches that employ a convex relaxation of the signal sparsity 
constraint. In X-ray computed tomography (CT), we propose 
an automatic scheme for extracting the convex hull of the 
inspected object from the measured sinograms; the obtained 
convex hull is used to capture the object contour information. 
We compare the proposed mask reconstruction schemes with 
the existing large-scale sparse signal reconstruction methods via 
numerical simulations and demonstrate that, by exploiting both 
the geometric contour information of the underlying image and 
sparsity of its wavelet coefficients, we can reconstruct this image 
using a significantly smaller number of measurements than the 
existing methods. 

I. Introduction 

Compressive sampling exploits the fact that most natural 
signals are well described by only a few significant (in mag- 
nitude) coefficients in some [e.g. discrete wavelet transform 
(DWT)] domain, where the number of significant coefficients 
is much smaller than the signal size. Therefore, for an p x 1 
vector x representing the signal and an appropriate p x p 
sparsifying transform matrix we have x = & s, where 
s = [si, «2? • • • •> Sp] T is an p x 1 signal transform-coefficient 
vector with most elements having small magnitudes. The idea 
behind compressive sampling or compressed sensing is to 
sense the significant components of s using a small number 
of linear measurements: 

y=$x (1) 

where y is an N x 1 measurement vector and $ is a known 
N x p sampling matrix with N < p\ here, we focus on 
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the scenario where the measurements, signal coefficients, and 
sampling and sparsifying transform matrices are real-valued. 
Practical recovery algorithms, including convex relaxation, 
greedy pursuit, and probabilistic methods, have been proposed 
to find the sparse solution to the underdetermined system (Q]), 
see Q for a survey. 

Compressive sampling takes the advantage of the prior 
knowledge that most natural signals are sparse in some 
transform domain. In addition to the signal sparsity, we use 
geometric constraints to enhance the signal reconstruction 
performance. In particular, we assume that the contour of 
the object under inspection is known and that the signal 
outside the contour is zero. A convex relaxation method was 
outlined in [2] for image reconstruction with both sparsity and 
object contour information. (Note that does not provide 
sufficient information to replicate its results and, furthermore, 
the method's development in El eqs. (4)-(6)] clearly con- 
tains typos or errors.) Here, we propose (i) iterative hard 
thresholding and convex relaxation algorithms that incorporate 
the object's contour information into the signal reconstruction 
process and (ii) an automatic scheme for extracting the convex 
hull of the inspected object (which captures the object contour 
information) from the measured X-ray computed tomography 
(CT) sinograms. 

We introduce our measurement model in Section HI and the 
proposed iterative hard thresholding methods in Section [TTT1 
Our mask convex relaxation algorithms are described in Sec- 
tion [TV] The experimental results are given in Section [VlJ 

We introduce the notation: || • || p and " T " denote the £ p norm 
and transpose, respectively, and the sparse thresholding oper- 
ator %(s) keeps the r largest-magnitude elements of a vector 
s intact and sets the rest to zero, e.g. 75 ([0, 1, — 5, 0, 3, 0] T ) = 
[0, 0, — 5, 0, 3, 0] T . The largest singular value of a matrix H is 
denoted by pn and is also known as the spectral norm of H. 
Finally, I n and nx i denote the identity matrix of size n and 
the n x 1 vector of zeros, respectively. 

II. Measurement Model 

We incorporate the geometric constraints via the following 
signal model: the elements of the p x 1 signal vector x = 

[xi,x 2 , • • • ,x p ] T are 

f [9a]i, i <EM 
1 " \ 0, i £ M Uj 



for i = 1, 2, . . . ,p, where [W s]i denotes the ith element of 
the vector & s, the mask M is the set of pm < P indices 
corresponding to the signal elements inside the contour of 
the inspected object, s is the p x 1 sparse signal transform- 
coefficient vector, and \P is the known orthogonal sparsifying 
transform matrix satisfying 



V 1 V 



(3) 



Therefore, the pm x 1 vector of signal elements inside the 
mask M (xi, i G M) is Xm = ^m,: s , where the pm x P 
matrix &m,-. contains the pm rows of & that correspond to 
the signal indices within the mask M. If the resulting $m, : 
has zero columns, the elements of s corresponding to these 
columns are not identifiable and are known to be zero because 
they describe part of the image outside the mask M. Define the 
set of indices I of nonzero columns of \Pm,-. containing p\ <p 
elements and the corresponding pi x 1 vector si of identifiable 
signal transform coefficients under our signal model. Then, 



x M = &mi si 



(4) 



where the pm xpi matrix #m,i is the restriction of \Pm,-. to the 
index set I and consists of the p\ nonzero columns of #m,:- 
Now, the noiseless measurement equation (Q} becomes [see 
also © and ®] 



y = (Px = $. M #m,i*i 



(5) 



where the N x pm matrix <P :j m is the restriction of the full 
sampling matrix $ to the mask index set M and consists of the 
Pm columns of the full sampling matrix $ that correspond to 
the signal indices within M. We now employ © and formulate 
the following constrained residual squared error minimization 
problem that incorporates both the geometric information (i.e. 
the knowledge of the inspected object's contour) and the signal 
sparsity constraint: 



(P ) : mm\\y - H si\\l subject to ||si|| 



< r 



(6) 



where ||si||o counts the number of nonzero elements in the 
vector si and H = <P : ,m &m,i- We refer to r as the signal 
sparsity level and assume that it is known. Finding the exact 
solution to © involves a combinatorial search and is therefore 
intractable in practice. In the following, we present greedy 
iterative schemes that aim at solving ©. 

III. Mask IHT and Mask DORE 

We first introduce a mask iterative hard thresholding (mask 
IHT) method and then propose its double overrelaxation 
acceleration termed mask DORE. 

Assume that the signal transform coefficient estimate is 
available, where q denotes the iteration index. Iteration (# + 1) 
of our mask IHT scheme proceeds as follows: 



= T r {s[ q) + ^ H T (y -Hs[ q) )) 



(7) 



where > is a step size chosen to ensure monotonically 
decreasing residual squared error, see also Section IIII-A1 
Iterate until and do not differ significantly. Upon 



convergence of this iteration yielding construct an 

estimate of the signal vector xm inside the mask M using 
#m,i In [3], we consider © with constant fi^ (not 

a function of q) set to /i^ = l/p%. For the full mask 
M = {1,2, and constant © reduces to the 

standard iterative hard thresholding (IHT) algorithm in (4). 

We now propose our mask DORE iteration that applies two 
consecutive overrelaxation steps after one mask IHT step to 
accelerate the convergence of the mask IHT algorithm. These 
two overrelaxations use the identifiable signal coefficient esti- 
mates sj^ and from the two most recently completed 
mask DORE iterations. Iteration (q + 1) of our mask DORE 
scheme proceeds as follows: 

1. Mask IHT step. 

s, = = %(s[ q) + ^ H T (y -Hs[ q) )) (8) 

where ji^ > is a step size chosen to ensure monotonically 
decreasing residual squared error, see also Section IIII-A1 

2. First overrelaxation. Minimize the residual squared error 
1 1 y — H sj || 2 with respect to si lying on the straight line 
connecting si and s[ 9 ^: 

zi = sj + ai (Si - s[ q) ) 

which has a closed-form solution: 

(Hsi-Hs^fiy-Hs!) 



(9a) 



(9b) 



3. Second overrelaxation. Minimize the residual squared error 
\y — H si || I with respect to si lying on the straight line 
connecting z\ and s| g_1 ^: 

which has a closed-form solution: 

(H z, - H (V - H z,) 



(10a) 



Q-1 



WHz.-Hst^g 



(10b) 



4. Thresholding. Threshold z\ to the sparsity level r: s\ = 

5. Decision. If ^ - Hsig < \\y -Hsi\\l assign = 
si; otherwise, assign = sj and complete Iteration q+l. 

Iterate until and sj^ do not differ significantly. As 

before, upon convergence of this iteration yielding 
construct an estimate of the signal vector xm inside the mask 
Musing ^ M ,iSi (+00) . 
A. Step size selection 

In Iteration 1 of our mask DORE and mask IHT schemes, 
we seek the largest step size p^ that satisfies 

Wv-hmI^Wv-Hs^w* (11) 

where si = si(s[°\ fi^) is computed using ([5]) with q = 0. 
We achieve this goal approximately as follows: Start with an 
initial guess for p^ > 0, compute the corresponding si = 

si( Sl (0) ,/i(°)), and 



if (ITU holds for the initial step size guess, double 

(repeatedly, if needed) until the condition (TTTT) for 

the corresponding si = si(s[°\ /jl^) fails; 

shrink (repeatedly, if needed) /i^ ^ by multiplying it with 

0.9 until (fTTT) for the corresponding si = 5i(sj , /i^) 

holds; 

complete Iteration 1 by moving on to Steps 2-5 in mask 



DORE or setting s\ 



si in mask IHT. 



In each subsequent Iteration q + 1 (q > 0), start with /i 

compute the corresponding si = /i^) in ([5]), 

and 



if 



W l 



(12) 
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does not hold for the initial step size fi 

/.(<?) 



(9-1) 



shrink /i^J by multiplying it (repeatedly, if needed) with 
0.9 until (Tf2l) for the corresponding si = si(s[ q \ fi^) 
holds; 

• complete Iteration q + 1 by moving on to Steps 2-5 in 
mask DORE or setting = si in mask IHT. 

Therefore, our step size /i^ is a decreasing piecewise constant 
function of the iteration index q. The step size obtained 
upon convergence (i.e. as q /* +00) is larger than or equal to 
0.9 /p 2 H , which follows easily from Theorem [T] below. 
Theorem I: Assuming that 



< ^ < l/p 2 H 



(13) 



and that the signal coefficient estimate in the q-th iteration s\ 
belongs to the parameter space 



(?) 



S r = {8ieW*: H|o<r} 



(14) 



then (IT21) holds, where si = si(sj , fi^) in (Tf2l) is computed 
using ([5]). Consequently, under the above conditions, the mask 
IHT and mask DORE iterations yield convergent monotoni- 
cally nonincreasing squared residuals \\y — H s^Wl as the 
iteration index q goes to infinity. 

Proof: See the Appendix. □ 



IV. Mask Convex Relaxation Methods 

Consider a Lagrange-multiplier formulation of © with the 
£q norm replaced by the l\ norm: 



(Pi): 



min(± \\y - 
6 i 



ll«i||i) (15) 



where r is the regularization parameter that controls the signal 
sparsity; note that the convex problem (1131) can be solved in 
polynomial time. Here, we solve (TT3T) using the fixed-point 
continuation active set (FPCas) an d gradient-projection for 
sparse reconstruction with debiasing methods in and O, 
respectively. We refer to these methods as mask FPCas an d 
mask GPSR, respectively. 



Fig. 1. Geometry of the parallel-beam X-ray CT system. 

V. Automatic Mask Generation from X-ray CT 
Sinograms Using a Convex Hull of the Object 

In X-ray computed tomography (CT), accurate object con- 
tour information can be extracted automatically from the 
measured sinograms. In particular, we construct a convex hull 
of the inspected object by taking intersection of the supports 
of the projections (over all projection angles) in the spatial 
image domain. 

To illustrate the convex hull extraction procedure, consider 
a parallel-beam X-ray CT system. Denote the measured sino- 
gram by po(t), where is the projection angle and t is the 
distance from the rotation center O to the measurement point. 
To obtain sufficient data for reconstruction, the range of t 
must be sufficiently large so that both ends of every projection 
Po(t) are zero. Define the range of the sinogram at angle 
by [a<9, be] = inf {[a, b] e R : pe(t) = for all t £ [a, b]} and 
the corresponding range in the spatial image domain: 

Aq = { (x, y) G R 2 : x cos 6 + y sin G [clq, be] } 

We construct the convex hull of the inspected object by taking 
the intersection H^Lo ^e- In practice, only a finite number K 
of projections is available at angles 0i, 62, . . . , Ok G [0, 7r), 
and the corresponding convex hull of the object can be 
computed as H^Li ^e k - Clearly, the angles #i, 62, • . • , 0k 
determine the tightness of the obtained convex hull. 

When imaging objects whose mass density is relatively high 
compared with that of the air, it is easy to determine the 
supports of the projections from the measured sinograms and 
extract the corresponding convex hull. For low-density objects 
such as pieces of foam, we need to choose carefully a threshold 
for determining these supports. 

VI. Numerical Examples 

In the following examples, we use the standard filtered back- 
projection (FBP) method |7, Sec. 3.3], which ignores both the 
signal sparsity and geometric object contour information, to 
initialize all iterative signal reconstruction methods. The mask 



DORE and DORE methods employ the following convergence 
criteria: 

n p+1) - Sl ^\\l/ Pl <e, \\ s (P+V-sW\\l/p<e (16) 

respectively, where e > denotes the convergence threshold. 

Shepp-Logan phantom reconstruction. We simulated 
limited-angle parallel-beam projections of an analog Shepp- 
Logan phantom with 1° spacing between projections and miss- 
ing angle span of 25°. Each projection is computed from its 
analytical sinogram using (8| function ellipse_sino . m] 
and (7) and then sampled by a receiver array containing 511 
elements. We then compute FFT of each projection, yielding 
N = 512 frequency-domain measurements; the corresponding 
frequency-domain sampling pattern is shown in Fig. |2(a)| 

Fig. |2(b)| depicts both the full and outer- shell masks of 
the phantom that we use to implement the DORE, GPSR, 
FPCas and mask DORE, GPSR, and FPC A s methods, re- 
spectively. Because of the nature of X-ray CT measurements, 
our full mask has circular shape containing p = 205859 
signal elements. The elliptical outer- shell mask containing 
Pm = 130815 w 0.6355 p pixels has been constructed from the 
phantom's sinogram using Hfc^i (fc-i)/i80> see Section IVl 
this choice of the mask implies that we have prior information 
about the shape of the outer shell of the Shepp-Logan phan- 
tom beyond the information available from the limited-angle 
projections that we use for reconstruction, see Fig. |2(a)[ 

Our performance metric is the peak signal-to-noise ratio 
(PSNR) of a reconstructed image x = [xi, #2, • • • , x p ] T inside 
the mask M: 



PSNR (dB) = 10 log 



[(max iGM Xi) - (min iGM Xi)f 



10 



Xi) 2 /p_ 



'M 



where x is the true image. 

We select the inverse Haar (Daubechies-2) DWT matrix 
to be the orthogonal sparsifying transform matrix SP\ the 
true signal vector s consists of the Haar wavelet transform 
coefficients of the phantom and is sparse: 



7866 ^ 0.0382 p. 



For the above choices of the mask and sparsifying transform, 
the number of identifiable signal transform coefficients is p\ = 
132450 « 0.6434 p. Note that ||s|| = ||si||o < Pi, implying 
that the identifiable signal coefficients are sparse as well. 
We compare the reconstruction performances of 

• mask DORE (r = 7000) and DORE (r = 8000) with 
e = 10 -14 [see (H6b l, where r are tuned for good PSNR 
performance; 

• the mask FPC A s, mask GPSR, FPC A s, and GPSR 
schemes, all using the regularization parameter r = 
10 ~ 5 ||i^ T i/||oo tuned for good PSNR performance; 

• the standard FBP method. 

(Here, we employ the convergence threshold tolP = 10 ~ 5 
for the mask GPSR and GPSR schemes, see |6l.) 

Figs. |2(c)] - |2(i)| show the reconstructions of various methods. 
To facilitate comparison, we employ the common gray scale to 



represent the pixel values within the images in Figs. |2(c)f|2(i)[ 
Clearly, taking the object's contour into account improves the 
signal reconstruction performance. 

Industrial object reconstruction. We apply our proposed 
methods to reconstruct an industrial object from real fan-beam 
X-ray CT projections. First, we performed the standard fan- 
to-parallel beam conversion (see [7, Sec. 3.4]) and generated 
parallel-beam projections with 1° spacing and measurement 
array size of 1023 elements, yielding N = 1024 frequency- 
domain measurements per projection. Our full mask has circu- 
lar shape containing p = 823519 signal elements. The outer- 
shell mask containing pm = 529079 « 0.6425 p pixels has 
been constructed from the phantom's parallel-beam sinogram 



using n fc= i A 



see Section N\ 



The mxm orthonormal sparsifying matrix \P is constructed 
using the inverse Daubechies-6 DWT matrix. 

We consider two measurement scenarios: no missing angles, 
i.e. all 180 projections available, and limited-angle projections 
with missing angle span of 20°, i.e. 160 projections available. 

We compare the reconstruction performances of mask 
DORE (r = 15000) and DORE (r = 20000) with e = 10" 8 ; 
the mask FPCas an d FPCas schemes using the regularization 
parameter r = 10 -6 \\H T y\\oo\ me standard FBP method. 
The reconstructions of mask FPCas an d FPCas are ver Y 
similar to those of mask DORE and DORE; hence we present 
only the mask DORE and DORE reconstructions in this 
example. Figs. |3(a)f|3(c)l show the reconstructions of the 
FBP, DORE, and mask DORE methods from 180 projections 
whereas Figs. |3(d)f|3(f)| show the corresponding reconstruc- 
tions from 160 limited-angle projections. Figs. |3(g)] - |3(i)| show 
the corresponding reconstruction profiles for slices depicted in 
Figs. |3(a)f|3(f)[ Observe the aliasing correction and denoising 
achieved by the sparse reconstruction methods. 

Appendix 

We now prove Theorem [T] Consider the inequality: 

||y - Hs^ HI - ||y - H^Wl = \\y - Hs^ f 2 - \\y - Hs^l 



> ||y - if 3k||! + -L||5i - s^\\l - \\H (5k - Sl M)\\l 
-\\y-Hs4l (Ala) 



(«)M|2 



1 



>(-^-P 2 h)\\si-s^\\1 



(Alb) 

where (lAlab follows by using the fact s\ in (EJi minimizes 

^||y-ff*i||i+||«i-«i^||l-A* (,) ||ff(«i-*i (a) )||l (A2) 

over all si G S r , see also (fT4l) . To see this, observe that (IA2b 
can be written as 

\\ Sl - s[ q) - H T (y - Hs[ q) )\\ 2 2 + const (A3) 

where const denotes terms that are not functions of sj. Finally, 
(lAlbl) follows by using the Rayleigh-quotient property O 
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Fig. 2. (a) 155 limited-angle projections in the 2-D frequency plane, (b) the full and outer-shell masks of the Shepp-Logan phantom, (c) FBP (PSNR = 
19.9 dB), (d) DORE (PSNR = 22.7 dB), (e) GPSR (PSNR = 22.9 dB), (f) FPC AS (PSNR = 22.5 dB), (g) mask DORE (PSNR = 25.8 dB), (h) mask 
GPSR (PSNR = 25.3 dB), and (i) mask FPC AS (PSNR = 26.4 dB) reconstructions. 



\Si - Si 



(9)1 



< 



Theorem 21.5.6]: \\H (si — si K 
Therefore, in each iteration, \\y — H s\ q) \\2 is guaranteed to 
not increase if the condition (ITSt holds. Since the sequence 
\\y — Hs^Wl is monotonically non-increasing and lower 
bounded by zero, it converges to a limit. 
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Fig. 3. FBP, DORE, and mask DORE reconstructions from (a)-(c) 180 projections and (d)-(f) 160 limited-angle projections; (g)-(i) the corresponding FBP, 
DORE, and mask DORE reconstruction profiles for slices depicted in (a)-(f). 
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